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HIGHLY  RESOLVED  SUB-TERAHERTZ  VIBRATIONAL  SPECTROSCOPY  OF 
BIOLOGICAL  MACROMOLECULES  AND  BACTERIA  CELLS 


1.  INTRODUCTION 

Sub-terahertz  (sub-THz)  spectroscopy  has  become  recognized  as  a  fast  technique 
for  sensing  and  detecting  biological  molecules  and  microorganisms.  It  is  based  on  the  interaction 
of  terahertz  (THz)  radiation  with  internal  molecular  vibrations  of  low  energy.  In  the  literature,  a 
majority  of  THz  experimental  data  have  been  reported  on  frequencies  above  1  THz  and  for 
relatively  small  biological  molecules  that  are  often  prepared  in  crystalline  form  (1 ).  Some 
progress  (past  4-6  years)  has  been  reported  in  the  literature  for  detecting  large  biological 
molecules  and  microorganisms  using  sub-THz  (below  1  THz)  spectroscopy  with  improved 
sensitivity  (2).  Molecular  dynamic  (MD)  simulations  for  relatively  small  molecular  components 
correlate  well  with  the  experimental  data  (2). 

Bacteria  are  very  complex  biological  objects.  Because  of  their  small  size  and 
relatively  low  absorption  coefficient,  THz  radiation  propagates  through  the  entire  object,  which 
allows  the  genetic  material  and  proteins  to  contribute  to  the  THz  signature  of  individual  bacteria 
or  spores.  In  the  past,  transmission  vibrational-mode  intensities  of  1-2%  could  be  obtained  when 
solid  samples  were  characterized  using  a  Fourier  transform  infrared  (FTIR)  spectrometer  (IFS 
66v;  Bruker  Daltonics;  Bremen,  Germany)  with  a  moderate  spectral  resolution  of  0.25  cm”1  (2). 
Raw  transmission  spectra  and  absorption  coefficient  data  were  calculated  for  the  quantitative 
characterization  of  biological  materials;  however,  discriminative  capabilities  needed  to  be 
improved.  Observed  spectroscopic  features  from  biological  cells  and  spores  are  caused  by  the 
interactions  between  THz  radiation  and  biological  macromolecules.  These  features  can  include 
combinations  of  low-energy  vibrational  modes  or  groups  of  modes  at  close  frequencies 
(vibrational  bands)  within  the  molecular  components  of  bacterial  cells  and  spores,  with  a 
significant  contribution  from  the  DNA  (2,  3).  The  entire  mechanism  of  MDs  is  still  not 
understood  and  needs  to  be  studied.  The  suggested  range  of  MD  relaxation  times  (x)  for 
processes  without  biological  molecular  conformational  change  varies  from  approximately  1.5  to 
650  ps  in  different  studies  ( 4 ,  5). 

The  width  of  individual  spectral  lines  and  the  intensity  of  resonance  features 
observed  in  sub-THz  spectroscopy  are  sensitive  to  the  relaxation  processes  of  atomic  dynamics 
(displacements)  within  a  macromolecule.  Those  relaxation  processes  determine  the 
discriminative  capabilities  of  sub-THz  spectroscopy.  Appropriate  spectral  resolution  must  be 
used  in  THz  spectroscopy  to  be  able  to  acquire  qualitative  as  well  as  quantitative  genetic 
information  used  to  identify  the  characteristic  signatures  of  biological  objects  that  will,  in  turn, 
increase  detection  accuracy  and  selectivity.  The  narrow  spectral  lines  are  due  to  the  existence  of 
long-lasting  dynamic  processes.  This  was  confirmed  by  the  relaxation  dynamics  of  side  chains  in 
macromolecules  that  were  observed  using  time-resolved  fluorescence  experiments  ( 6 ),  which 
required  the  application  of  highly  resolved  vibrational  THz  spectroscopy.  The  majority  of 
available  THz  spectroscopy  systems  do  not  have  the  high  spectral  resolution  that  is  required  to 
observe  narrow-width  spectral  lines. 


1 


Experimental  measurements  with  good  spectral  resolutions  have  been 
demonstrated  (7,  8)  in  which  intense  and  narrow  spectral  features  were  observed  from  biological 
molecules  and  bacteria  with  spectral  widths  between  0.05  and  0.2  cm-1.  These  features  were  not 
evident  in  previous  experimental  measurements  that  used  a  spectral  resolution  of  0.25  cm-1. 
These  findings  suggested  the  coexistence  of  diverse  relaxation  dynamics  mechanisms  that  are 
relevant  to  the  sub-THz  frequency  region  (9).  Multiple  intense  and  specific  resonance  signatures 
may  provide  a  way  for  the  development  of  reliable  discriminative  capabilities  using  the  new 
sub-THz  frequency-domain  spectroscopy  described  herein.  This  was  not  possible  using  the 
methods  discussed  in  the  previous  literature  studies.  The  controversies  found  in  the  relatively 
young  sub-THz  spectroscopy  field  were  due  to  experimental  problems,  such  as  poor  spectral 
resolution  and  sample  preparation  techniques  that  were  not  adequate  for  the  problem. 

The  Vibratess  (Charlottesville,  VA)  sub-THz  spectrometer  is  a  continuous-wave, 
frequency-domain  spectroscopic  THz  instrument  with  imaging  capability.  The  detector  operates 
at  room  temperature  without  the  need  for  cryogenic  cooling.  Frequency  transmission  spectra 
were  obtained  in  the  sub-THz  region  between  315  and  480  GHz  for  macromolecules  and 
biological  species.  Due  to  the  high  sensitivity,  good  spectral  resolution,  and  a  spatial  resolution, 
the  spectroscopic  instrument  allows  for  the  observation  of  intense  and  narrow  spectral 
resonances  in  transmission  and  absorption  spectra  for  nanogram  amounts  of  samples  obtained 
from  biological  materials  with  spectral  line  widths  as  narrow  as  ~0.1  cm-1.  Experimental  spectra 
from  the  new  sub-THz  spectrometer,  obtained  using  nanogram  quantities  of  biological 
macromolecules  and  biological  cells  and  spores,  were  compared  to  the  MD  simulations. 

The  majority  of  THz  systems  that  are  commercially  available  use  the  time-domain 
approach.  Some  disadvantages  of  this  approach  are:  (1)  fast  Fourier  transform  (FFT)  is  a 
powerful  technique,  but  there  are  limitations  on  the  use  of  FFT  with  discrete  data  points  (not 
continuous  function);  (2)  multiple  frequencies  radiate  onto  the  sample,  which  leads  to  unreliable 
absorptions  and  measurements  (e.g.,  cross-vibrations);  and  (3)  collection  of  limited  data  points 
will  force  averaging  of  absorption  peaks. 

The  frequency-domain  THz  system  used  in  this  study  has  the  following 
advantages:  (1)  new  and  reliable  THz  signature,  (2)  better  resolution  with  frequency-domain 
THz  signatures,  and  (3)  simulation  that  can  be  used  to  easily  explain  any  single  THz  peak  of 
interest. 


2.  EXPERIMENTAL  SETUP 

2.1  New  Sub-THz  System 

The  new  Vibratess  sub-THz  system  (VsTHz;  Figure  1)  has  several  advantages 
over  the  commercial  systems  available  at  the  time  of  this  study  (e.g.,  Bruker  IFS  66v).  The 
VsTHz  system  does  not  require  milligram  quantities  of  biological  material,  a  liquid  helium 
cooled  detector,  operating  under  vacuum,  or  purging  with  dry  gas.  Although  some  specificity 
was  observed  in  different  bacterial  spectra  using  the  Bruker  IFS  66v  instrument,  the  absorption 
peak  intensities  for  Bacillus  subtilis  and  Escherichia  coli  were  rather  close.  This  suggests  that  the 


2 


spectral  resolution  of  0.25  cm  '1  did  not  provide  a  sufficient  level  of  discrimination.  Better 
spectral  resolution  is  required  for  better  discrimination  and  clear  specificity.  Other  important 
characteristics  of  a  sub-THz  system  are  small  sample  size  requirements,  a  simple  sample- 
preparation  procedure,  room-temperature  operation,  good  reproducibility,  and  user-friendly 
operation  for  wide-spread  adoption  of  a  new  system. 

The  new  VsTHz  system  was  developed  with  improved  spectral  resolution,  which 
was  close  to  1  GHz  (0.03  cm-1).  THz  systems  described  in  previous  literature  lacked  spectral 
resolution.  The  VsTHz  system  operates  in  a  frequency  domain  at  room  temperature  and  with 
local  enhancement  of  the  electromagnetic  field;  therefore,  coupling  of  the  THz  radiation  with  the 
sample  biomaterials  was  increased  (8, 10-14). 

The  VsTHz  system  operates  in  a  frequency  range  between  315  and  480  GHz,  with 
an  average  output  power  of  ~0.2  mW.  It  includes  a  sub-THz  electronically  tunable  source  that  is 
based  on  Schottky  frequency  multipliers  (Virginia  Diodes,  Inc.  [VDI];  Charlottesville,  VA)  and  a 
submicron,  precision  motorized  stage  to  provide  three-dimensional  (3D)  scanning  capability  for 
a  detector  subsystem,  relative  to  a  6  mm  wide,  multichannel  sample  chip  holder  (Figure  2).  The 
detector  system  has  a  microdetector  (Schottky  diode  from  VDI)  with  a  replaceable  beam- lead 
microprobe  (Lichtenberger  Consulting;  Charlottesville,  VA)  and  a  planar  microdetector  circuit 
(75),  which  is  mounted  in  a  custom  microwave  guide  housing. 


Figure  1.  Vibratess  spectroscopic  sensor  prototype.  Micro-THz  detector,  microchannels  sample 
chip  holder,  and  objective  for  optical  visualization  (horizontal  camera). 
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Figure  2.  Microchannels  (along  the  y  axis  from  -3000  to  3000  |am)  sample  chip  (6  mm  wide). 


The  sample  chip  holder  (Figure  2)  is  constructed  of  an  array  of  microchannels  in  a 
gold  or  copper  layer  on  a  substrate.  The  substrate  is  any  material  with  a  low  absorption 
coefficient  and  a  low  refractive  index  that  is  in  the  THz  range.  This  configuration  generates  the 
edge  effect,  which  is  necessary  for  enhancement  of  the  electromagnetic  field  of  the  THz 
radiation  and  provides  regions  for  holding  biological  materials.  The  substrate  material  we  used 
was  polymethylmethacrylate  with  channel  widths  of  5-20  pm.  The  thickness  of  the  gold  layer 
(2-5  pm)  determines  the  channel  depth.  A  high-magnification  optical-imaging  system  with  a 
long  working-distance  objective  allows  accurate  positioning  of  components.  The  waveguide  and 
microdetector  circuit  were  designed  using  a  high-frequency  structure  simulator  to  separate  THz 
radiation  from  a  direct  current  Schottky  diode  signal.  Spatial  resolution  is  better  than  200  pm  and 
is  restricted  by  the  opening  size  of  the  commercially  available  microdetector  waveguide  (150  x 
200  pm2).  In  the  future,  smaller  microdetector  waveguides  can  be  used  to  increase  the  spatial 
resolution  (lower  than  200  pm). 

2.2  Experimental  Procedure 

Care  was  taken  to  prevent  the  system  from  overheating  and  to  ensure  the 
temperature  stabilization  of  the  microchannel  chip  and  protein  samples  under  investigation.  We 
took  steps  to  minimize  the  heat  from  the  illumination  source  and  ensured  that  the  temperature  in 
the  room  was  controlled  at  24-25.5  °C  during  measurements. 

The  minimum  possible  amount  of  sample  material  in  a  channel  used  for 
interrogation  was  0.05  ng  in  an  area  of  10  x  200  pm,  with  a  channel  depth  of  3  pm.  A  protein 
solution  droplet  size  of  0.2  pL,  with  concentration  of  0.25  mg/mL,  yielded  enough  solution  to 
spread  on  an  area  ~1  x  2  mm2.  This  was  more  than  the  amount  required  to  provide  high 
sensitivity  for  THz  characterization.  It  is  crucial  to  ensure  that  the  sample  material  layer  is  very 
thin.  In  our  example,  the  dried  layer  thickness  during  the  sample  spectrum  scan  was  ~2.5  pm, 
which  was  slightly  less  than  the  channel  depth.  A  thicker  sample  layer  (>3  pm)  containing  water 
will  have  a  significantly  higher  dielectric  constant  in  our  frequency  range  as  compared  with  air 
and  plastic.  This  will  modify  the  optical  path  for  radiation  through  the  sample  chip  as  compared 
with  the  background  of  empty  channels,  which  will  result  in  artificial  values  of  the  calculated 
transmission  above  1.  In  addition,  a  very  thin  substrate  (only  25  pm  thick)  bends  under  the 
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weight  of  the  sample  material  solution,  which  also  causes  poor  stability  and  reproducibility.  In 
our  case,  the  reproducible  results  were  obtained  after  ~2  h  of  material  drying.  Additional 
improvements  are  needed  to  reduce  this  waiting  time  from  hours  to  several  minutes  after  sample 
solution  deposition  on  the  chip. 

2.2.1  Sample  Preparation  and  Deposit 

For  this  study,  a  micropipette  was  used  to  deliver  a  0. 1-0.3  pL  droplet  of 
biological  sample  in  liquid  suspension  (80%  water  and  20%  methanol)  on  a  preselected  spot  on 
the  sample  chip  (array  of  microchannels)  in  the  sample  holder.  Within  a  minute,  the  droplet  on 
top  of  the  sample  chip  shrank  quickly  in  size  and  then  evaporated  completely.  The  microdetector 
probe  was  positioned  several  microns  above  the  droplet’s  spot  on  the  sample  chip.  Only  ~20  ng 
of  biomaterial  is  required  for  the  VsTHz  system,  as  compared  with  the  milligram  sample  size 
required  for  other  THz  systems  that  are  mentioned  in  the  literature. 

2.2.2  Data  Collection  and  Measurements 

Measurements  and  data  collection  were  performed  10  min  after  the  droplet  was 
deposited  and  the  sample  was  dried.  The  acquired  data  consisted  of  an  array  containing  the 
frequency-dependent  signal  and  other  pertinent  metadata  including  data,  time,  and  the  3D  probe 
coordinates.  A  spectral  resolution  >1  GHz  was  achieved,  along  with  high  sensitivity.  A  typical 
signal-to-noise  ratio  obtained  was  — 103  (depending  on  frequency)  at  room  temperature  (without 
the  requirement  for  air  evacuation  or  purging  with  dry  nitrogen). 

2.2.3  Characterization  of  the  Sub-THz  System 

Figure  3  shows  the  periodicity  of  the  transmitted  signal  at  a  frequency  of 
15.5  cm"1  as  the  probe  was  scanned  across  the  sample  chip  microchannels  from  a  blank  run  (no 
biological  sample  deposited).  As  expected,  it  showed  the  same  periodicity  as  the  array  itself, 
although  the  data  was  slightly  offset  in  the  x  direction  from  the  channels. 
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Scanning  across  the  channels  (YIG  Voltage  =4.1  V) 


Figure  3.  Periodicity  of  a  THz  signal  (measured  signal  in  microvolts,  gray  line) 
in  the  direction  ( x  in  micrometers  from  the  center  of  the  sample  chip  [Figure  2]) 
perpendicular  to  the  microchannels  (gray  bar). 


Figure  4  shows  the  transmitted  signal  at  a  frequency  of  15.5  cm”1  of  a  sample  of 
B.  subtilis  at  four  microchannel  y-axis  values  (-250,  -300,  -350,  -400  pm  [Figure  2])  as  the 
probe  was  scanned  across  the  sample  chip  microchannels  (in  the  x-axis  direction).  Figure  2 
shows  the  microchannel  x  and  y  directions.  The  sample  was  located  between  -450  and  150  pm 
on  the  microchannel  x  axis,  but  a  longer  scan  was  used,  and  the  observed  transmission  was  close 
to  1  as  the  probe  was  scanned  across  empty  microchannels  (no  biological  sample).  Transmission 
was  reduced  to  about  60%  in  area  as  the  probe  was  scanned  across  microchannels  with  biological 
sample.  The  actual  reduction  depended  on  the  amount  of  biological  material  in  each  channel. 
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X  coordinate,  jim 

Figure  4.  Transmission  scan  across  channels  with  B.  subtilis  cell  material. 


3.  THZ  MOLECULAR  DYNAMIC  (MD)  SIMULATION  OF  BIOLOGICAL 

MATERIALS 

In  a  parallel  study,  a  computational  modeling  technique  was  developed  to 
simulate  vibrational  spectra  of  proteins  using  the  energy  minimization,  normal  mode  analysis, 
and  MD  approaches.  The  MD  simulations  were  performed  using  the  assisted  model  building 
with  energy  refinement  (AMBER)  simulation  program  package  (University  of  California;  San 
Francisco,  CA)  as  described  by  Alijabbari  et  al.  (3),  with  a  modified  procedure  to  improve  the 
convergence  of  modeling  results  for  proteins. 

3.1  Insulin  MD  Simulation 

To  model  the  insulin  molecule,  it  was  necessary  to  correctly  choose  an  insulin 
crystal  structure  file  from  the  protein  data  bank  (PDB)  that  most  closely  corresponded  to  an 
experimental  material  used  in  our  absorption  measurements.  The  PDB  database  contains  436 
structures  related  to  human  insulin.  The  experimental  material  (Sigma- Aldrich  Corporation 
[St.  Louis,  MO],  catalog  no.  10908  SIGMA)  consisted  of  a  mix  of  insulin  and  Zn2+  ions 
(Figure  5A-C).  It  is  known  that  three  insulin  molecules  are  able  to  form  a  stable  complex  with 
Zn2+  ion  in  the  center  (1)  as  shown  in  Figure  5A  (16). 

The  PDB  file  (3W7Z)  was  used  for  MD  simulation  of  the  entire  complex.  The 
crystalline  structure  of  this  complex  is  shown  in  Figure  5A.  Other  structures  of  insulin  can  be 
simpler  and  contain  only  one  insulin  molecule  that  consists  of  two  chains  with  disulfide  bridges, 
as  shown  in  Figure  5B,C.  Two  different  PDB  files  were  used  as  starting  points  to  simulate  a 
single  insulin  molecule — one  of  them  was  PDB  file  3W7Z,  which  uses  Zn2+  ions  to  stabilize  the 
structure,  and  the  second  was  PDB  file  3I3Z.  This  file  (3I3Z),  which  represents  only  one 
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molecule,  was  chosen  because  it  contains  the  atomic  coordinates  of  human  insulin  crystallized 
without  Zn2+  ions. 

As  shown  in  Figure  5A,  three  insulin  molecules  formed  a  stable  complex  with  the 
Zn2+  ion  (green)  in  the  center  through  interaction  between  Zn2+,  water  molecules,  and  histidine 
residues  (data  from  3W7Z  PDB  file).  In  Figure  5B,  two  chains  of  insulin  molecules  with  one 
intra-  and  two  intermolecular  disulfide  bridges  (yellow).  Insulin  was  crystallized  with  Zn2+  ions 
(data  from  3W7Z  PDB  file).  In  Figure  5C,  in  contrast  with  the  molecule  in  Figure  5B,  the 
molecule  was  crystallized  without  Zn2+  ions  (data  from  3I3Z  PDB  file).  Atom  positions  were 
slightly  different  but  similar  in  configuration. 


A 


Three  insulin  molecules 
form  a  stable  complex 
with  the  Zn2+  ion  (green) 
in  the  center.  The  complex 
is  formed  through 
interaction  between  Zn2+, 
water  molecules,  and 
histidine  residues  (data 
from  3W7Z  PDB  file). 


B 

Two  chains  of  insulin 
molecules  with  one  intra- 
and  two  intermolecular 
disulfide  bridges  (yellow). 
Insulin  was  crystallized  with 
Zn2+  ions  (data  from  3W7Z 
PDB  file). 


/ 


C 

In  contrast  with  a  molecule 
in  Figure  5B,  this  molecule 
was  crystallized  without 
Zn2+  ions  (data  from  3I3Z 
PDB  file).  Atom  positions 
are  slightly  different  but 
similar  in  configuration. 


Figure  5.  Three  insulin  molecules. 


3.2  Lysozyme  MD  Simulation 

There  are  several  3D  molecular  structures  for  lysozyme  available  in  the  PDB, 
including  1BWH,  1LYS,  and  4QEQ  PDB  structures  determined  from  X-ray  analysis.  It  is 
important  to  choose  the  best  structure  that  most  closely  correlates  with  the  real  structure  of 
purchased  lysozyme  for  MD  simulations  and  analysis.  In  previous  literature  studies,  the  structure 
1BWH  was  found  to  have  a  low  resolution.  The  4QEQ  PDB  structure  has  been  obtained  only 
recently.  The  quality  of  X-ray  crystal  structures  of  4QEQ  and  1LYS  PDB  can  be  found  at  the 
PDB  website  for  lysozyme  files  (17, 18). 

The  information  provided  in  this  section  was  based  on  the  global  validation 
metrics  for  the  4QEQ  and  1LYS  lysozyme  PDB  files.  According  to  this  global  validation  metric 
criteria,  the  4QEQ  structure  has  better  Ramachandran  (3)  and  Sidechain  outlier  scores,  which 
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means  that  the  structure  of  4QEQ  is  of  better  quality.  Several  MD  simulations  were  conducted 
using  the  1BWH  PDB  file.  This  included  three  different  initial  random  atomic  velocities  and 
total  energy  (Etot)  of  a  molecular  system  at  room  temperature  (293  K)  and  at  an  elevated 
temperature  of  298  K. 

4.  DISCUSSION 

The  ability  to  quickly  and  reliably  discriminate  between  different  bacterial 
species,  toxins,  and  proteins  using  sub-THz  spectroscopy  would  provide  significant  benefits  to 
the  medical  first  responders  and  soldiers.  For  example,  in  the  medical  field,  it  would  enable  a 
faster  and  more-tailored  treatment  for  a  patient  or  affected  area  once  a  bacterial  organism  or 
toxin  is  identified  as  the  cause  of  an  infection,  spillage,  release,  and  so  on. 

To  demonstrate  the  capabilities  of  the  VsTHz  spectrometer,  and  to  fully 
understand  the  basic  fundamentals  of  THz  absorption  by  biomaterials,  the  study  started  with 
commonly  used  proteins  that  had  peptide  deformylase  structures  in  the  public  database  before 
beginning  a  THz  study  of  bacterial  cells  and  some  of  their  molecular  components  (DNA  and 
thioredoxin).  Extensive  experimental  characterization  and  MD  simulations  of  THz  spectroscopic 
signatures  were  conducted  using  human  insulin  protein  and  chicken  egg-white  lysozyme  protein. 

Proteins  with  different  structures  will  absorb  and  vibrate  in  response  to  light  of 
different  frequencies.  A  protein  can  be  considered  like  a  tree  with  different  size  leaves  and 
branches;  as  air  passes  over  the  tree,  the  leaves  and  branches  vibrate  and  slightly  rotate  at 
different  frequencies.  Protein  end  groups  are  like  the  leaves  and  amino  acids  are  like  the 
branches  of  the  tree. 

4.1  Insulin  THz  Study 

4.1.1  Insulin  THz  Experimental  Spectra 

The  human  insulin  protein  (10908,  molecular  weight  -5800  Da)  was  purchased 
from  Sigma-Aldrich  (http://www.sigmaaldrich.com/catalog/product/sigma/i0908?lang= 
en&region=US).  Human  insulin  is  a  two-chain  polypeptide  hormone  produced  by  the  P-cells  of 
pancreatic  islets.  The  a  and  P  chains  are  joined  by  two  interchain  disulfide  bonds.  The  a  chain 
contains  an  intrachain  disulfide  bond.  Insulin  regulates  the  cellular  uptake,  utilization,  and 
storage  of  glucose,  amino  acids,  and  fatty  acids  and  inhibits  the  breakdown  of  glycogen,  protein, 
and  fat. 


The  human  insulin  protein  we  obtained  was  actually  a  mixture  of  insulin  and  Zn2+ 
ions.  Through  direct  communication  with  Sigma-Aldrich  technical  support,  it  was  found  that 
there  was  no  available  PDB  file  that  could  be  used  to  characterize  the  3D  structure  of  Sigma- 
Aldrich  insulin. 
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The  transmission  spectra  shown  in  Figure  6  was  obtained  by  following  the 
experimental  setup,  sample  preparation,  and  data  collection  procedures  discussed  in  Sections  2.2 
to  2.2.2.  Figure  6  shows  the  transmission  change  over  time  while  the  material  fully  dried  (80, 
130,  and  250  min).  In  this  example,  the  sample  material  was  diluted  to  very  low  concentration  of 
0.25  mg/mL  to  slow  the  drying  process.  The  data  show  that  the  position  of  all  transmission 
minima  on  the  frequency  scale  was  preserved,  although  the  absolute  value  of  transmission  was 
changing. 
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Figure  6.  Transmission  change  depending  on  drying  time  (insulin  concentration  of  0.25  mg/mL 

in  50%  ethanol  and  droplet  size  ~1  pL). 


Transmission  spectra  were  recalculated  to  absorbance,  A  =  -log  (7),  where  T  is 
transmission.  Figure  7  compares  the  absorbance  of  dried  sample  (transmission  at  250  min,  the 
red  line  plot  on  Figure  6)  with  the  averaged  results  from  measurements  at  different  spots  on  a 
chip.  Once  again,  the  frequencies  of  most  transmission  minima  were  reproduced.  The  red  spectra 
shown  in  Figures  6  and  7  are  the  same;  Figure  6  is  shown  in  transmission  scale,  and  Figure  7  is 
shown  in  absorbance  scale.  The  red  spectrum  in  Figure  7  was  chosen  over  the  average  spectrum 
of  all  spots  for  comparison  with  MD  simulations  because  it  was  more  reliable  and  had  stronger 
maximum  absorptions. 
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Figure  7.  Variability  of  experimental  spectra  that  depends  on  the  radiation  intensity  variation  in 
different  channels  and  spots  on  a  microchip.  (Insulin  concentration  of  0.25  mg/mL  in  50%  ethyl 

alcohol  and  droplet  size  of  0.3  pL). 


Figure  8  demonstrates  the  absorbance  spectra  from  thicker  sample  layers.  In  this 
case,  the  drying  process  required  several  days  to  be  completed.  The  results  were  less  reliable 
because  the  background  and  sample  measurements  were  conducted  on  different  days  to  allow 
sufficient  time  for  the  material  to  fully  dry. 
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Figure  8.  Measured  absorbance  spectrum  from  thick  insulin  sample  after  drying  for  several  days. 
(Insulin  concentration  of  0.5  mg/mL  and  droplet  size  of  0.3  pL). 
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As  shown  in  Figure  8,  the  most-important  result  was  that  the  dominant  absorption 
peak  was  identified  at  frequencies  15.5-15.8  cm"1.  The  most-reliable  data,  shown  by  the  red  lines 
in  Figures  6  and  7,  were  further  used  for  comparison  with  results  from  MD  simulations. 

4.1.2  Insulin  THz  MD  Simulations 

To  model  the  insulin  molecule,  it  was  necessary  to  choose  a  PDB  insulin  crystal 
structure  that  most  closely  corresponded  to  an  experimental  material  used  in  our  absorption 
measurements.  The  PDB  file  (3W7Z)  contained  information  about  the  complex  of  three  insulin 
molecules  with  Zn2+  ions;  therefore,  it  was  used  for  the  MD  simulation  of  the  entire  complex. 

The  crystalline  structure  of  this  complex  is  shown  in  Figure  5A.  The  material  in  the  experiment 
was  amorphous  and  did  not  have  the  crystalline  structure  shown  in  Figure  5A,  which  was 
specifically  built  for  X-ray  analysis.  The  unit  in  this  amorphous  material  can  be  simpler  and  may 
contain  only  one  insulin  molecule,  which  consists  of  two  chains  with  disulfide  bridges,  as  shown 
in  Figures  5B  and  5C.  Two  different  PDB  files  were  used  as  a  starting  point  to  simulate  a  single 
insulin  molecule:  PDB  file  3W7Z,  which  used  Zn2+  ions  to  stabilize  the  structure,  and  PDB  file 
3I3Z.  The  PDB  file  3I3Z,  which  represents  only  one  molecule,  was  chosen  because  it  contains 
atomic  coordinates  of  human  insulin  crystallized  without  Zn2+  ions. 

To  characterize  the  theoretical  spectrum  of  insulin,  and  to  better  analyze  and 
understand  the  experimental  spectra,  we  conducted  multiple  simulations  with  two  different  goals. 
The  first  goal  was  to  estimate  the  variability  of  the  spectra  due  to  the  difference  of  initial  atomic 
coordinates  taken  from  crystal  structures  by  modeling  one  insulin  molecule.  For  this  purpose,  we 
modeled  single  insulin  molecules  that  were  crystallized  with  Zn2+  ions  (3W7Z  PDB  file)  and 
without  Zn2+  (3I3Z  PDB  file). 

The  second  goal  was  to  estimate  the  effect  of  weak  intermolecular  bonds  between 
insulin  molecules  by  modeling  a  complex  of  three  insulin  molecules  together  with  Zn2+  ions. 

This  question  was  addressed  by  comparing  the  simulation  results  obtained  for  a  single  molecule 
with  those  of  the  complex  taken  from  the  same  PDB  file,  3W7Z. 

MD  simulations  were  performed,  and  absorption  coefficients  were  found  using 
the  procedures  described  in  Section  3.1.  A  detailed  description  of  the  Vibratess  protocol 
approach  for  MD  simulation  using  the  AMBER  program  can  be  found  in  a  study  done  by  Globus 

o 

et  al.  (2).  We  simulated  a  complex  containing  a  biomolecule  and  an  8  A  shell  of  TIP3P  water.  In 
the  preparation  of  this  complex,  a  constant  number,  volume,  and  temperature  (NVT)  ensemble 
was  used  to  raise  the  temperature  to  293  K. 

The  system  was  heated  for  ~16  ps,  and  the  protein  atoms  were  restrained  using  a 
10  kcal/mol/A"  force  constant.  During  this  heating  process,  the  bonds  involving  hydrogen  were 
fixed.  Constant  pressure  periodic  boundary  conditions  were  used  to  scale  the  system  volume 
during  100  ps  to  reach  a  density  of  ~1  g/cm3.  In  these  procedures,  a  10  A  real-space  cutoff  was 
used  with  a  2  fs  integration  time  step.  After  the  system  attained  the  selected  values  of 
temperature  and  density,  random  velocities  from  the  Maxwellian  distribution  were  assigned  to  all 
atoms,  followed  by  another  equilibration  step;  then  a  constant  number,  pressure,  and  temperature 
(NPT)  ensemble  was  used  for  further  energy  minimization.  After  equilibration,  a  600  ps 
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production  run  was  performed  in  a  constant  number,  volume,  and  energy  (NVE)  ensemble  to 
avoid  undesirable  effects  on  atomic  motions  because  of  coupling  to  an  external  thermal  bath. 
During  the  simulation,  the  coordinates  of  all  atoms  of  the  system  were  recorded  every  20  fs  for 
the  entire  production  run.  Atomic  trajectories  collected  in  MD  simulations  were  converted  to  the 
covariance  matrix  of  atomic  random  displacements,  <RiRk>.  The  force-field  matrix  is  found  in  a 
quasi-harmonic  approximation  using  the  relation  between  the  covariance  matrix  and  the  inverse 
of  the  force-constant  matrix  (<RiRk>  =  kisT  [F  -  l]ik),  where  R  is  the  random  displacement  value 
between  atom  i  and  atom  k,  T  is  temperature,  and  F  is  a  force  constant  (3,  4).  Diagonalization  of 
the  F  matrix  gives  eigenfrequencies  (normal  mode  frequencies)  and  eigenvectors  (displacement 
vectors-normal  modes). 

Using  atomic  trajectories  from  the  constant  energy  and  volume  MD  simulations, 
insulin  sub-THz  vibrational  spectra  and  absorption  coefficients  were  calculated  in  a  quasi¬ 
harmonic  approximation.  The  absorption  coefficient  spectra  a(v),  as  a  function  of  the  frequency 
(V),  can  be  calculated  through  the  relationship  between  a  and  the  imaginary  part  of  dielectric 
permittivity  (5) 


a(v)  =  Wv2^ 


k 


Sk 

v2)2+W2v2 


(l) 


where  Vk  is  the  normal  mode  frequency,  calculated  by  diagonalization  of  the  force-constant 
matrix,  and  5k  is  the  oscillator  strengths  computed  for  all  vibrational  modes,  k.  The  spectral  line 
width,  W  =  0.12  cm-1,  which  is  close  to  the  width  of  spectral  features  observed  in  experiments 
with  the  Vibratess  spectrometer,  was  used  for  all  of  our  vibrational-mode  simulations  in  sub-THz 
range. 


Better  simulation  convergence  and  improved  consistency  between  simulated 
vibrational  frequencies  and  experimental  data  were  obtained  using  a  new  procedure  for 
averaging  mass- weighted  covariance  matrices  of  atomic  trajectories  in  MD  simulations  (2). 
Averaging  of  only  six  matrices  gave  consistent  results.  The  overall  stability  of  protein  MDs  was 
estimated  using  a  comparison  of  two  spectra  that  were  calculated  with  a  covariance  matrix  from 
an  entire  production  run  and  from  an  averaged  covariance  matrix  ( 3 ,  20-22). 

Figures  9-12  show  the  MD  simulation  results  compared  with  the  experimental 
spectrum  (red  curve  in  Figure  7).  The  MD  simulations  also  confirmed  that  THz  absorption  peaks 
have  a  peak  width  of  approximately  0.2  cm”1.  This  means  that  the  proper  THz  system  should 
have  a  THz  resolution  of  at  least  0.01  cm”1  (20  points  per  THz  absorption  peak  or  higher).  Better 
resolution  requires  a  longer  time  for  collecting  data  and  does  not  provide  any  additional 
information. 


Figure  9  shows  the  experimental  THz  spectrum  of  insulin,  compared  with  the  MD 
simulation  THz  spectrum  obtained  for  a  single  insulin  molecule  using  the  3W7Z  PDB  file.  There 
was  good  correlation  with  the  most-prominent  experimental  THz  absorption  peak 
(15.5-16.2  cm”1)  when  one  insulin  molecule,  crystallized  with  Zn2+,  was  modeled.  The  vertical 
green  lines  in  Figure  9  represent  the  THz  absorption  peaks  found  in  the  experimental  and 
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modeling  spectra  (six  matching  peaks).  The  most-prominent  MD  simulation  peak  is  located  at 
-15.9  cm-1.  The  most-prominent  experimental  THz  spectrum  peak  is  located  at  -15.8  cm"1. 


11  12  13  14  15  16 


Frequency,  cm' 

Figure  9.  Insulin  experimental  THz  spectrum  compared  with  the  MD  simulation  THz  spectrum 
for  a  single  insulin  molecule  using  the  3W7Z  PDB  file  with  Zn2+. 


Figure  10  shows  the  experimental  THz  spectrum  compared  with  the  MD 
simulation  THz  spectrum  using  a  single  insulin  molecule  from  the  3I3Z  PDB  file.  The  MD 
simulation  THz  spectrum  that  was  based  on  the  3I3Z  PDB  structure  crystallized  without  Zn2+; 
however,  the  most-prominent  peak  is  located  at  -14.6  cm1,  which  is  not  in  agreement  with  the 
experimental  THz  spectrum.  The  vertical  blue  lines  in  Figure  10  represent  the  THz  absorption 
peaks  found  in  experimental  and  modeling  spectra  (five  matching  peaks).  The  most-prominent 
experimental  THz  spectrum  peak  is  located  at  -15.8  cm"1. 
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Figure  10.  Experimental  THz  spectrum  compared  with  the  MD  simulation  THz  spectrum  for  a 
single  insulin  molecule  from  the  3I3Z  PDB  file  without  Zn2+. 


Figure  1 1  shows  the  experimental  THz  spectrum  compared  with  the  MD 
simulation  THz  spectrum  obtained  using  the  complex  of  three  insulin  molecules  and  a  Zn2+  ion. 
The  MD  THz  spectrum  simulated  using  the  3W7Z  PDB  structure  for  a  complex  of  three 
molecules  and  a  Zn2+  ion  (Figure  5A)  showed  good  correlation  with  experimental  absorption 
peaks  at  frequencies  11.4,  12.2,  14.2,  16.3,  and  16.6  cm”1,  as  shown  in  Figure  11.  The  vertical 
green  lines  in  Figure  1 1  represent  the  THz  absorption  peaks  found  in  experimental  and  modeling 
spectra  (eight  matching  peaks).  There  are  three  most-prominent  MD  simulation  THz  peaks 
located  at  -11.2,  13.2,  and  14.2  cm”1.  The  most-prominent  experimental  THz  spectrum  peak  is 
located  at  -15.8  cm”1. 
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Figure  1 1 .  Experimental  spectrum  compared  with  the  MD  simulation  THz  spectrum  for  the 
complex  of  three  insulin  molecules  and  Zn2+  ion  PDB  file. 


Figures  9  and  1 1  show  results  with  better  fits  to  the  experimental  THz  spectra  as 
compared  with  Figure  10.  Based  on  the  correlations  (vertical  lines)  shown  in  Figures  9  and  11,  it 
can  be  assumed  that  the  experimental  material  might  be  a  mixture  of  two  components:  (1)  a 
single  insulin  molecule  and  (2)  a  complex  of  three  insulin  molecules  crystallized  with  Zn2+.  A 
combined  MD  simulation  THz  spectrum  was  calculated  from  the  two  3W7Z  PDB  structure  MD 
simulations  where  two  spectra  were  taken  with  1:1  ratios.  This  new,  combined  spectrum 
(Figure  12)  had  an  improved  correspondence  as  compared  with  the  experimental  data. 

Figure  12  shows  the  experimental  THz  spectrum  as  compared  with  the  average 
MD  simulation.  The  THz  spectra  were  obtained  from  a  single  insulin  molecule  MD  simulation 
and  from  the  complex  of  three  insulin  molecules  MD  simulation  crystallized  with  a  Zn2+  ion 
(data  3W7Z  PDB  file).  The  vertical  green  lines  in  Figure  12  represent  the  THz  absorption  peaks 
found  in  the  experimental  and  modeling  spectra  (nine  matching  peaks). 
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Figure  12.  Insulin  experimental  THz  spectrum  compared  with  that  of  the  average  MD 
simulation.  The  vertical  green  lines  represent  the  THz  absorption  peaks  found  in  experimental 

and  modeling  spectra  (nine  matching  peaks). 


Results  from  the  experimental  THz  spectrum  and  those  from  various  MD 
simulations  on  human  insulin  in  the  sub-THz  frequency  range  indicate  that  the  origin  of  the 
vibrational  modes  is  heterogeneous;  in  other  words,  it  is  a  combination  of  multiple  structures  and 
ions. 

4.2  Lysozyme  THz  Study 

4.2.1  Lysozyme  THz  Experimental  Spectra 

Lysozyme  is  a  small,  globular  protein  composed  of  129  amino  acid  residues. 
Lysozyme  is  able  to  damage  bacterial  cell  walls  by  catalyzing  the  hydrolysis  of  peptidoglycans 
(found  in  the  cell  walls  of  bacteria,  especially  Gram-positive  bacteria)  and  hydrolyzing  the 
glycosidic  bond  that  connects  M acetyl muramic  acid  with  the  fourth  carbon  atom  of 
iV- acetyl  glucosamine.  Lysozyme  is  widely  used  in  the  enzymatic  lysis  of  microbial  cells  (23,  24). 

The  lysozyme  material  that  was  used  for  experimental  characterization  was 
purchased  from  Sigma- Aldrich  and  contained  L6876  lyophilized  powder  and  protein  >90% 
(>40,000  units/mg  of  protein).  The  lysozyme  sample  with  a  concentration  of  5  mg/mL  in 
distilled  water  was  prepared  at  the  U.S.  Army  Edgewood  Chemical  Biological  Center. 
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The  procedure  that  was  used  for  experimental  characterization  was  basically  the 
same  as  that  used  for  insulin  and  presented  in  Section  3.1.  The  previously  used  conventional 
light  source  that  caused  biosample  material  heating  during  measurements  was  substituted  with 
the  light-emitting  diode  base  source. 

Material  with  an  initial  concentration  of  5  mg/mL  was  diluted  to  a  concentration 
of  0.83  mg/mL  with  distilled  water  and  mixed  1:1  with  ethyl  alcohol.  Samples  with  the  initial 
and  diluted  concentrations  were  both  characterized. 

Figure  13  shows  the  background  THz  spectral  measurements  at  three  different 
locations  on  the  microchannels  sample  chip  (as  depicted  in  Figure  2)  but  very  close  together  (x- 
axis  locations  were  -448.55,  -449.30,  and  -449.98  pm  ).  These  results  showed  good 
reproducibility  between  the  three  background  THz  spectral  measurements  (red,  blue,  and  green 
lines)  and  the  average  spectrum  (black)  before  the  lysozyme  sample  was  deposited. 


> 

— j 


CU 

c 

O) 

CO 


Frequency,  cm-1 

Figure  13.  Reproducibility  and  variability  of  independent,  background  spectral  measurements  in 
the  same  position  (sample  chip  microchannel  [Figure  2]  x-axis  locations  were  -448.55,  -449.30, 
and  -449.98  pm)  before  the  lysozyme  sample  was  deposited. 


Figure  14  shows  the  comparison  between  the  average  background  THz  spectrum 
(black  line  in  Figure  13)  and  the  chicken  egg-white  lysozyme  THz  spectrum  (blue  line)  produced 
after  the  lysozyme  sample  was  deposited.  The  signal  intensity  difference  between  the  black  and 
blue  lines  represents  the  absolute  absorption. 
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Figure  14.  Comparison  between  average  background  spectrum  and  chicken  egg-white  lysozyme 

protein  spectrum. 


Figure  15  shows  the  experimental  variability  of  the  lysozyme  transmission  spectra 
at  three  very  close  locations  on  the  microchannels  sample  chip  (Figure  2);  x-axis  locations  were 
-448.18,  -449.53,  and  -449.75  pm,  which  are  shown  as  blue,  red,  and  green  lines,  respectively. 
The  black  THz  spectrum  represents  the  average  of  the  three  lysozyme  transmission  spectra.  The 
three  lysozyme  transmission  spectra  were  produced  using  the  average  background  THz  spectrum 
(black  line  in  Figure  13).  Note  that  the  experimental  THz  spectral  features  (positions  of  local 
maximums  and  minimums)  on  the  x  axis  (frequency)  are  reproducible,  at  better  than  0.02  cm-1. 
The  reproducibility  of  the  spectral  features  positions  (local  maximums  and  minimums)  on  the 
frequency  x  axis  (Figure  15)  was  good  at  >0.02  cm-1. 
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Figure  15.  Variability  of  lysozyme  transmission  spectra  at  three  very  close  locations  on  the 
microchannels  sample  chip  (Figure  2);  x-axis  locations  were  -448.18,  -449.53,  and  -449.75  pm. 
The  black  THz  spectrum  represents  the  average  of  the  three  lysozyme 

transmission  spectra. 


Figure  16  shows  the  change  in  the  THz  absorbance  spectra  based  on  the  amount 
of  sample  material  that  was  deposited  on  the  same  spot  in  the  sample  chip  microchannel. 
Frequencies  of  all  maxima  and  minima  were  reproduced  extremely  accurately,  but  the  absolute 
values  were  different  and  did  not  follow  Beer’s  law.  This  can  be  understood  because  the  amount 
of  sample  material  analyzed  was  restricted  by  the  height  of  the  sample  chip  holder  microchannel 
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Figure  16.  Lysozyme  THz  absorbance  spectra  for  different  amounts  of  material.  Frequencies  of 
all  maxima  and  minima  are  reproduced  accurately. 
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Figure  17  shows  that  the  larger  concentrations  (blue  and  brown  lines,  5  and 
0.41  mg/mL,  respectively)  did  not  reduce  the  transmission  on  the  same  location  (y  =  0  pm). 
However,  the  transmission  spectra  had  slight  variation  at  different  locations  on  the 
microchannels  (y  =  0  and  -150  pm,  blue  and  brown,  respectively).  Spectral  variations  are  due  to 
the  variations  in  the  sample  and  how  the  sample  spread  over  different  microchannels.  THz 
measurements  at  different  locations  or  spots  on  the  microchannels  chip  play  a  big  role  in  the 
measured  THz  spectra,  as  shown  by  the  blue  spectrum  at  a  location  y  =  0  pm  on  the  y  axis  as 
compared  with  the  brown  spectrum  at  a  location  y  =  -150  pm  on  the  y  axis.  Future  studies  will 
examine  the  role  of  sample  spread  and  amount  on  the  measured  THz  spectra  to  find  the  optimal 
experimental  settings.  The  transmission  spectra  in  Figure  17  need  to  be  converted  to  absorbance 
spectrum  (Figure  18)  to  be  further  compared  with  the  THz  MD  simulations.  Figure  18  shows  the 
average  absorbance  spectrum  of  lysozyme,  using  transmission  spectra  from  Figure  17. 


Figure  17.  Transmission  spectra  of  lysozyme  at  two  different  concentrations  (0.41  and  5  mg/mL, 
which  are  shown  as  green  and  blue  lines,  respectively)  and  at  two  different  y-axis  locations  on 
the  microchannels  sample  chip  (y  =  0  and  -150  pm,  which  are  shown  as  blue  and  brown  lines, 

respectively). 
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Figure  18.  Averaged  absorbance  spectrum  of  lysozyme  using  transmission  spectra 

from  Figure  17. 


4.2.2  Lysozyme  THz  MD  Simulation 

Because  the  lysozyme  molecule  was  much  larger  in  size,  as  compared  with 
insulin,  it  was  expected  to  produce  many  more  spectral  peaks  than  would  insulin,  and  these 
spectral  peaks  could  overlap. 

There  were  several  3D  molecular  structures  for  lysozyme  available  in  the  PDB 
file.  These  included:  (1)  1BWH,  (2)  1LYS,  and  (3)  4QEQ  PDB  structures,  which  were 
determined  by  X-ray  analysis.  It  was  important  to  choose  the  structure  that  most  closely 
correlated  with  that  of  the  purchased  and  analyzed  lysozyme  material.  The  structure  1BWH  had 
already  been  obtained  with  a  low-resolution  X-ray  system.  The  4QEQ  PDB  structure  was 
obtained  with  a  higher-resolution  X-ray  system.  The  X-ray  crystal  structures  of  4QEQ  and  1LYS 
can  be  found  at  the  PDB  website  (17,  18).  According  to  these  global  validation  metrics  criteria, 
the  4QEQ  structure  had  better  Ramachandran  (19)  and  sidechain  outlier  scores  than  did  1BWH 
or  1LYS,  which  meant  that  the  4QEQ  structure  was  of  better  quality. 

Several  MD  simulations  were  conducted  using  the  1BWH  PDB  structure  with 
three  different  initial  random  atomic  velocities  and  with  Etot  at  room  temperature  (293  K)  and  at 
an  elevated  temperature  of  298  K. 

For  example,  Figures  19  and  20  show  progress  in  the  simulation  predictive 
capability  for  lysozyme.  Figure  19  shows  the  lysozyme  MD  simulation  results  at  room 
temperature,  T  =  293  K.  Three  different  plots,  calculated  at  three  different  total  energies  of  a 
protein  (Etot  =  -50947,  -50915,  and  -50888  kcal/mol,  which  are  shown  as  blue,  pink,  and  green 
lines,  respectively,  in  Figure  19),  show  that  these  spectra  were  sensitive  to  the  Etot  values  used 
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even  though  experiments  were  carried  out  at  room  temperature.  This  could  be  expected  because 
the  lysozyme  melting  temperature  under  pH  7  is  only  71  °C  (25),  and  it  is  not  far  enough  from 
room  temperature.  This  means  that  a  change  in  the  measurement  temperature  can  strongly  affect 
the  spectrum.  For  this  reason,  we  also  conducted  a  simulation  not  only  at  T  =  293  K,  but  also  at 
T  =  298  K,  which  was  the  temperature  near  the  spectrometer  at  the  time  of  measurement. 

Figure  20  shows  the  lysozyme  absorption  spectrum  that  was  simulated  for  298  K. 


Figure  19.  Lysozyme  MD  simulation  results  at  room  temperature,  T  =  293  K.  Three  different 
plots,  calculated  at  three  different  total  energies  of  a  protein  (Etot  =  -50947,  -50915,  and 
-50888  kcal/mol,  which  are  shown  as  blue,  pink,  and  green  lines,  respectively). 


Four  different  MD  simulations  were  performed  at  different  Etot  values,  as  shown 
in  Figures  19  and  20.  The  MD  simulation  at  Etot  =  -52260  kcal/mol  shows  a  better  correlation 
with  the  experimental  measurement  (Figure  20)  than  with  the  other  three  spectra  (Figure  19). 
Although  the  correlation  between  MD  simulation  at  Etot  =  -52260  kcal/mol  and  the 
experimental  measurement,  as  shown  in  Figure  20,  showed  that  some  measured  absorption  peaks 
occurred  at  the  same  frequency  as  in  simulation,  the  results  were  not  favorable.  The  next  MD 
simulation  was  made  using  different  Etot  values  to  increase  correlation  results  and  to  get  close 
to  a  100%  match  with  all  experimental  peaks. 


23 


Freq u en cy ,  cm'1 

Figure  20.  Lysozyme  THz  MD  simulation  absorption  spectrum,  which  was  simulated  with  the 
PDB  file  1BWH  using  a  room  temperature  value  of  T  =  298  K  and  at  the  single  total  energy  of  a 

protein,  Etot  =  -50260  kcal/mol. 


In  addition,  as  was  the  case  with  insulin,  it  was  necessary  to  explore  more  than 
one  crystal  structure  from  those  available  in  the  PDB.  MD  simulation  modeling  was  conducted 
with  two  more  crystal  structures  using  two  new  PDB  files:  (1)  1LYS  and  (2)  4QEQ,  which 
provide  slightly  different  atomic  positions.  Figure  21  compares  the  two  THz  MD  simulations  that 
used  the  two  PDB  files  (4QEQ  and  1LYS,  which  are  shown  as  brown  and  pink  lines, 
respectively)  at  T  =  298  K  (without  Cl“  ions)  with  experimental  THz  spectra  (green  lines). 
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Figure  21.  Two  THz  MD  simulations  using  two  PDB  files  (4QEQ  and  1LYS,  which  are  shown 
as  brown  and  pink  lines,  respectively)  at  T  =  298  K  (without  Cl  ions)  with  experimental  TFIz 

spectra  (green  line). 


Although  many  peaks  in  the  simulated  and  experimental  data  have  the  same  or 
very  similar  frequencies,  MD  simulations  were  also  conducted  with  two  lysozyme  structures  and 
Cl“  ions  added  to  the  model  (system  of  protein  and  water)  because  these  ions  were  present  in  the 
experimental  material.  The  results  are  shown  in  Figures  22  and  23.  Figure  22  shows  the  THz 
absorption  spectra  from  an  MD  simulation  using  two  model  structures  (4QEQ  and  1LYS,  which 
are  shown  as  brown  and  pink  lines,  respectively)  with  CF  ions. 
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Figure  22.  THz  absorption  spectra  from  an  MD  simulation  using  two  model  structures  (4QEQ 
and  1LYS,  which  are  shown  as  brown  and  pink  lines,  respectively)  with  Cl-  ions. 


Figure  23  shows  a  comparison  between  two  THz  absorption  spectra,  one  from  an 
MD  simulation  using  PDB  file  4QEQ  with  Cl“  ions  (brown  line)  and  the  other  from  experimental 
THz  measurements  (T  =  298  K)  (green  line).  In  Figure  23,  the  lysozyme  MD  simulation  spectra 
used  an  averaging  correlation  matrix  interval  of  100  ps. 


Frequency,  cm  1 

Figure  23.  Comparison  between  two  THz  absorption  spectra  from  an  MD  simulation  using  PDB 
file  4QEQ  with  Cl-  ions  (dark  red  line)  and  from  experimental  THz  measurements  at  T  =  298  K 

(green  line). 
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As  expected  and  shown  Figure  23,  the  THz  absorption  spectrum  using  PDB  file 
4QEQ  with  Cl  ions  (brown  line)  shows  a  closer  correlation  with  the  experimentally  measured 
THz  spectrum  (green  line).  THz  MD  simulated  spectra  were,  however,  dependent  on  several 
parameters  that  were  not  yet  optimized.  Additional  efforts  sought  to  improve  the  convergence 
between  MD  simulation  and  experimental  spectra  by  using  optimizing  simulation  parameters,  as 
shown  in  Figure  24. 


Freq uency,  cm 

Figure  24.  Lysozyme  experimental  THz  spectrum  (green  plot)  and  THz  MD  spectrum  (blue  line) 
determined  using  averaging  correlation  matrices  for  five  intervals:  20,  30,  40,  50,  and  60  ps. 


Figure  24  shows  the  lysozyme  experimental  THz  spectrum  (green  line)  and  THz 
MD  spectrum  (blue  line)  that  were  determined  using  averaging  correlation  matrices  for  five 
intervals:  20,  30,  40,  50,  and  60  ps,  and  not  just  one  time  interval  (100  ps),  as  is  shown  in 
Figure  23  (brown  line).  Figure  24  shows  better  correlation  between  the  experimental  and  MD 
simulation  spectra  of  lysozyme.  The  two  red  oval  regions  indicate  the  main  discrepancies 
observed  in  this  evaluation.  Simulated  lysozyme  spectra  are  sensitive  to  the  initial  atomic 
position  provided  by  PDB  files.  It  was  important  to  explore  several  PDB  crystal  structures  and 
conduct  simulations  with  Cl“  ions  to  determine  whether  the  interaction  with  ions  was  responsible 
for  some  features  in  experimental  spectra  and  to  resolve  the  discrepancies.  It  was  found  that  the 
correlation  between  the  experimental  data  and  MD  simulation  resulted  in  peak  position 
differences  on  a  frequency  scale  of  <0.1  cm"-1  and  with  good  relative  peak  intensities. 

We  evaluated  the  dominant  THz  peak  at  about  13  cm-1,  shown  in  Figure  24,  to 
understand  the  primary  causes  of  vibrations,  as  detailed  in  Figure  25.  Figure  25a  shows  the  same 
graph  as  Figure  24,  with  two  highlighted  regions  of  discrepancies  (red  ovals)  between  the  MD 
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simulation  and  the  experimental  spectrum.  In  Figure  25b,  a  quasi-harmonic  analysis  of  the 
13  cm  1  resonant  frequency  reveals  that  arginine  128  and  glutamine  121  residues  have  the 
highest  oscillation  amplitudes  in  the  z  direction  as  compared  with  the  other  residues.  The  atoms 
that  make  the  biggest  contribution  to  the  oscillations  are  marked  by  spheres  (white  for  hydrogen, 
blue  for  nitrogen,  red  for  oxygen,  and  grey  for  carbon).  Figure  25c  shows  the  variation  over  time 
in  the  distance  between  glutamine  121  and  arginine  128  residues  in  the  z  direction.  The 

o 

amplitude  was  extremely  large  at  ~1  A. 


13.1  cm-1  frequency  in 
lysozyme  THz  spectrum 
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Figure  25.  Diagram  of  a  potential  method  to  speed  up  MD  simulations:  (a)  same  graphic  as  that 
of  Figure  24,  with  attention  on  the  dominant  vibrational  peak  at  about  13  cm-1;  (b)  quasi¬ 
harmonic  analysis  of  the  13  cm'1  resonant  frequency  reveals  that  arginine  128  and  glutamine  121 
residues  have  highest  oscillation;  and  (c)  variation  over  time  of  the  distance  between 
glutamine  121  and  arginine  128  residues  in  the  z  direction. 


MD  simulation  using  multiple  time  intervals,  as  shown  in  Figure  24,  improved  its 
correlation  with  the  experimental  THz  spectra  of  lysozyme.  It  also  provided  an  in-depth 
understanding  of  the  responsible  vibrations.  Unfortunately,  each  MD  simulation  takes  a  few  days 
to  generate.  Figure  25  shows  a  new  approach  that  can  be  used  to  speed  up  the  MD  simulation. 
Using  this  method,  individual  experimental  THz  spectroscopic  THz  vibrational  peaks  are 
interrogated  and  modeled,  but  the  entire  frequency  range  is  not.  In  Figure  25,  evaluation  of  the 
13  cm  1  resonant  frequency  revealed  that  arginine  128  and  glutamine  121  residues  had  the 
highest  oscillation  amplitudes  in  the  z  direction  as  compared  with  the  other  residues.  Arginine 
128  and  glutamine  121  were  the  primary  causes  of  the  13  cm-1  peak,  regardless  of  the  chosen 
PDB  file  structure  (4QEQ  or  1LYS).  Fine  tuning  the  13  cm-1  peak  was  accomplished  much  more 
rapidly  by  using  MD  simulation  for  only  the  arginine  128  and  glutamine  121  with  different 
settings,  rather  than  by  using  the  previous  MD  simulation  approach  in  which  all  atoms  and 
entities  in  the  structure  were  analyzed.  This  new  approach  will  significantly  reduce  the  MD 
simulation  time  and  improve  our  understanding  of  the  underlying  mechanisms  that  produce 
experimental  THz  peaks. 
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4.3 


Bacteria  THz  Study 


4.3.1  Bacteria  THz  Experimental  Spectra 

A  microsyringe  was  used  to  deposit  the  B.  subtilis  var.  niger  spore  suspension  in 
the  correct  location  and  amount.  Unfortunately,  when  the  0.05  pL  droplet  of  bacteria  in 
suspension  dried,  it  formed  a  “coffee  ring”  as  shown  in  Figure  26.  Most  of  the  bacteria 
suspended  in  the  droplet  moved  to  the  edge,  which  left  a  ring  of  bacteria  after  drying.  Therefore, 
the  THz  detector  tip  had  to  be  moved  to  the  edge  of  the  dried  droplet  for  this  study,  which 
required  the  use  of  a  second  vertical  camera  to  provide  a  top  view  for  precise  ring-edge  location. 


Electron  micrograph  of 
0.05ul  droplet  of  2mg/ml 
subtilis  spores  on  THz 
Size  740  x  430  micrometers 


Figure  26.  Drying  behavior  for  a  0.05  pL  droplet  of  2  mg/mL  of  B.  subtilis  spores  on  a  THz 
microchannel  sample  chip  holder  (evaporated  to  leave  a  ring  similar  to  a  “coffee-ring”  effect). 


The  E.  coli  cell  studies  presented  in  the  next  several  paragraphs  are  found  in  an 
article  by  Globus  et  al.  and  are  used  to  show  the  capability  of  the  VsTHz  system  with  bacterial 
cells  (9). 


To  demonstrate  the  capabilities  of  the  VsTHz  spectrometer,  transmission  spectra 
from  bacterial  cells  and  some  of  their  molecular  components  (DNA  and  thioredoxin)  were 
measured  and  recalculated  for  the  absorption  coefficient  per  mass  of  material  as  described  by 
Globus  et  al.  (2).  The  reliability  of  the  transmission  results  was  verified  by  measuring  samples 
with  different  amounts  of  material  and  at  different  coordinates  in  the  microchannel  array. 

Figure  27  shows  the  transmission  spectra  of  E.  coli  cells  in  two  different  amounts,  1  and  2  pL. 
The  two  spectra  show  good  correlation  (reproducibility)  between  the  frequency  locations  of  the 
spectral  features  and  the  transmission.  In  addition,  the  absorption  was  generally  higher  (lower 
transmission  values)  for  the  2  pL  sample,  which  indicated  that  the  absorption  was  from  the 
sample. 
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Figure  27.  Transmission  of  E.  coli  cells  using  two  samples  with  different  amounts  of  material. 


Figure  28  shows  the  absorbance  spectra  from  one  droplet  (blue  line)  and  two 
droplets  (brown  line)  of  B.  subtilis  spores.  Each  droplet  size  was  0.05  pL,  contained  0.2  pg  of 
B.  subtilis  material,  and  spread  in  an  area  of  0.3  mm  in  diameter.  The  B.  subtilis  sample 
concentration  was  4  mg/mL. 


Figure  28.  Absorbance  spectra  of  B.  subtilis  spores  from  0.05  pL  (blue  line,  1  droplet)  and 
0.1  pL  (brown  line,  2  droplets).  B.  subtilis  spore  sample  concentration  was  4  mg/mL. 
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Transmission  and  absorption  spectra  of  E.  coli  DNA  are  shown  in  Figures  29  and 
30.  Very  narrow  spectral  lines  were  detected,  with  widths  of  ~0. 1-0.2  cm-1.  The  peak  positions 
were  preserved  in  the  spectra  of  samples  having  different  amounts  of  material,  but  the  intensities 
of  the  absorption  peaks  calculated  per  unit  mass  were  often  reduced  for  larger  sample  amounts 
(Figure  30). 


Figure  29.  Transmission  spectrum  of  E.  coli  DNA  with  500  ng  of  material  in  the  droplet. 


We  assumed  that  the  thickness  of  the  sample  material  after  drying  (at  the  edge  of 
the  ring)  was  proportional  to  its  mass,  which  was  determined  from  the  volume  of  solution  and 
the  concentration  of  material  in  solution.  This  gave  us  the  mass  of  solid  material  in  the  drop  after 
evaporation  of  water.  We  further  assumed  that  the  spot  area  was  the  same  when  more  material 
was  added.  However,  at  the  time  of  this  study,  we  could  not  control  the  spot  size  and  could  not 
be  sure  that  this  area  was  the  same.  This  may  explain  the  apparent  dependence  of  absorption 
coefficient  on  mass,  as  shown  in  Figure  30. 
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Figure  30.  Absorption  coefficient  spectra  of  E.  coli  DNA.  Scaling  was  not  observed  with  the 

amount  of  sample  material  studied. 


As  shown  in  this  and  in  earlier  work  by  Globus  et  al.  (2),  spectral  lines  are  sharper 
in  thin  samples,  and  they  have  significant  broadening  and  damping  of  spectral  features  in  bulk 
samples  that  is  probably  because  of  less-oriented  material.  Very  thick  samples  usually  do  not 
provide  well-resolved  spectral  features  in  the  sub-THz  region.  This  effect  can  explain  the 
sensitivity  of  sub-THz  spectroscopic  characterization  to  the  specific  sample  preparation 
techniques. 


Further  development  of  the  sample  delivery  approach  is  expected  to  provide 
increased  reproducibility  and  efficiency  because  the  sample  presentation  in  the  spectrometer 
detection  region  will  be  more  accurate  and  controlled. 

Figure  31  compares  the  averaged  transmission  spectra  of  E.  coli  cells,  as 
measured  on  the  Bruker  FTIR  spectrometer  (IFS  66v),  with  the  results  from  the  VsTHz 
spectrometer.  Different  transmission  scales  were  used  for  the  two  instruments,  and  a  significantly 
reduced  amount  of  sample  was  used.  Overlapping  spectral  peaks  (narrow  and  intense  spectral 
lines)  in  the  spectrum,  measured  by  the  Bruker  FTIR  instrument,  are  resolved  by  the  VsTHz 
spectrometer. 
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Figure  31.  Transmission  spectra  of  E.  coli  cells  measured  with  two  different  instruments,  the 
Bruker  FTIR  (IFS  66v,  red  line)  and  the  VsTHz  frequency-domain  spectroscopic  sensor 

(blue  line). 


There  was  a  shift  in  the  major  transmission  minimum  position  from  14.6  to 
14.9  cm"1.  This  0.3  cm"1  shift  to  a  higher  frequency  was  explained  by  the  different  orientations 
of  the  cells  in  the  two  setups,  which  were  needed  because  spectra  are  sensitive  to  the  orientation 
of  biomaterial  with  respect  to  the  electric  field  vector  in  THz  radiation  (2). 

Comparison  of  the  parameters  and  conditions  used  with  the  new  VsTHz 
spectrometer  with  those  of  the  Bruker  spectrometer  IFS  66v  in  the  sub-THz  operational  range 
(11-17  cm"1)  has  demonstrated  an  order  of  magnitude  better  spectral  resolution  and  10-15  times 
higher  peak  intensities. 

The  high  sensitivity  of  the  VsTHz  spectrometer  allowed  2  orders  of  magnitude 
less  sample  to  be  deposited  on  the  sample  chip  (yielding  a  spot  of  1-2  mm2).  The  VsTHz 
detector  waveguide  interrogates  and  acquires  signals  on  only  a  fraction  of  a  1-2  mm2  sample 
spot  because  the  detector  antenna  is  only  0.12  x  0.15  mm  in  size. 

Because  overlapping  of  individual  resonances  was  still  possible,  the  estimates  for 
the  dissipation  factor  from  the  width  of  spectral  lines  in  Figures  28-30  provided  an  upper  limit 
for  y  between  0.1  and  0.2  cm"1. 
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4.3.2 


Bacteria  THz  MD  Simulation 


Unfortunately,  because  any  bacterial  cell  consists  of  many  biomolecules  and 
chemicals,  and  because  PDB  files  do  not  exist  for  bacteria,  an  MD  simulation  could  not  be 
performed.  But  part  of  a  cell’s  biography  can  be  modeled  as  described  in  this  section. 

Figure  32  shows  the  experimental  spectra  of  the  E.  coli  protein  thioredoxin,  with 
computational  modeling  results  using  MD  simulations  and  a  damping  factor  of  y  =  0.12  cm-1. 
These  results  were  used  to  confirm  the  existence  of  the  observed  narrow  and  intense  resonance 
features  in  the  sub-THz  transmission  and  the  absorption  spectra  of  biological  materials,  as 
measured  with  VsTHz  spectrometer. 
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Figure  32.  Absorption  spectra  of  protein  thioredoxin  from  E.  coli ,  obtained  using  MD  simulation 
and  experimental  results  from  the  VsTHz  spectrometer. 


The  computational  modeling  procedure  was  conducted  exactly  as  it  was  described 
in  detail  by  Alijabbari  et  al.  (3).  The  only  difference  was  that  a  smaller  value  of  the  dissipation 
factor  (0.12  cm-1)  was  used  in  this  study,  which  was  close  to  the  approximate  experimental 
widths  of  spectral  features  observed  with  the  VsTHz  spectrometer.  Not  all  peaks  (as  shown  in 
Figure  32)  were  reproduced  in  the  measured  and  simulated  spectra.  Simulation  parameters  have 
not  yet  been  optimized,  and  the  same  value  of  y  was  used  to  calculate  absorption  for  all 
vibrational  modes. 
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Initial  MD  simulation  and  experimental  THz  spectra  confirmed  the  existence  of 
intense  and  narrow  absorption  lines,  which  can  be  used  for  discrimination  between  different 
bacteria  and  strains.  The  narrow  absorption  lines  were  confirmed  in  an  article  by  Qiu  et  al.  ( 6 ) 
from  x  scales  in  thioredoxin  from  femtosecond-resolved  fluorescence  spectroscopy  experiments. 
The  aforementioned  authors  ( 6 )  also  observed  longer  quenching  dynamic  processes  with  time 
scales  of  275-615  ps  at  a  hydrogen  bond  distance,  which  can  give  local  fluctuations  with 
vibration  spectral  line  widths  of  <0.1  cm-1. 


5.  CONCLUSIONS 

The  VsTHz  spectrometer  (frequency-domain  THz  spectroscopic  system 
prototype,  operating  at  room  temperature)  was  tested  with  a  variety  of  biological  samples.  This 
prototype  system  provided  high  spectral  resolution  and  significantly  improved  detection 
sensitivity  and  reliability  in  the  sub-THz  operational  range  as  compared  with  commercially 
available  spectrometers.  The  results  of  this  study  demonstrated  highly  resolved  transmission 
(absorption)  spectra  from  only  10-20  ng  of  biological  macromolecules,  bacterial  cells,  and 
spores. 


The  experimental  results  revealed  very  intense  and  narrow  spectral  features  from 
biological  molecules  and  bacteria  with  widths  of  ~0. 1-0.2  cm1.  The  results  provided  completely 
new  information  about  the  interactions  between  THz  radiation  and  biological  materials,  which 
confirmed  the  diverse  relaxation  dynamic  mechanisms  that  are  relevant  to  sub-THz 
spectroscopy.  This  demonstration  of  multiple  intense  and  specific  resonance  features  provides 
potential  conditions  for  the  reliable  discriminative  capabilities  of  this  technology  to  the  level  of 
bacterial  strains  in  extremely  small  sample  volumes.  However,  further  evaluation  is  required  to 
fully  assess  the  limits  of  these  capabilities. 

Experimental  results  and  improvements  on  the  existing  VsTHz  spectrometer  will 
create  a  foundation  for  the  development  of  new  types  of  advanced  biological  THz  systems.  In 
addition,  this  new  instrument  can  be  used  for  monitoring  interactions  between  biomaterials  and 
reagents  and  for  studying  conformational  changes  and  biomedical  processes  in  near  real  time. 

Improvements  that  are  important  for  improved  computational  modeling  include 
the  following:  finding  the  PDB  structures  that  best  represent  the  experimental  material;  more 
accurately  representing  components  of  the  protein  surroundings  such  as  water,  ion  molecules, 
and  so  on;  and  averaging  the  results  from  several  runs  while  varying  the  equilibration 
parameters.  Finally,  using  quasi-harmonic  analysis  of  individual  resonant  frequencies,  and  using 
that  information  to  optimize  the  MD  simulation  can  result  in  greatly  improved  simulation  times. 

Future  experimental  efforts  should  include  further  investigation  regarding  the  role 
of  sample  thickness,  location,  water  content,  and  freshness  (changing  the  absorption  spectrum  of 
the  sample  material  due  to  degradation).  In  addition,  future  efforts  should  also  strive  to  increase 
the  number  of  experimental  runs  on  each  sample  to  ensure  a  good  statistical  representation  of  the 
in-class  diversity.  Finally,  statistical  methods  such  as  “Principal  Component  Analysis”  should  be 
used  to  quantify  the  correlations  between  the  MD-simulated  and  experimental  THz  spectra. 
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